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ABSTRACT 

We study the role of dust-dust collisional charging in protoplanetary discs. We 
show that dust-dust collisional charging becomes an important process in determining 
the charge state of dust and gas, if there is dust enhancement and/or dust is fluffy, so 
that dust surface area per disc volume is locally increased. 

We solve the charge equilibrium equations for various disc environments and dust 
number density 77, using general purpose graphic processors (GPGPU) and CUDA 
programming language. We found that as dust number density 77 increases, the charge 
distribution experience four phases. In one of these phases the electrostatic field E 
caused by dust motion increases as i? cx ry^. As a result, macroscopic electric discharge 
takes place, for example at 77 = 70 (in units of minimum-mass solar nebula (MMSN) 
values, considering two groups of fluffy dust with radii 10~^ cm, 10^ cm). We present 
a model that describes the charge exchange processes in the discs as an electric circuit. 
We derive analytical formulae of critical dust number density for lightning, as functions 
of dust parameters. 

We estimate the total energy, intensity and event ratio of such discharges ('light- 
ning'). We discuss the possibility of observing lightning and sprite discharges in proto- 
planetary discs by Astronomically Low Frequency [ALF) waves, IR images, C/T^ lines, 
and high energy gamma rays. We also discuss the effects of lightning on chondrule 
heating, planetesimal growth and magnetorotational instability of the disc. 

Key words: methodsmumerical — planetary systems:formation — planetary sys- 
tems:protoplanetary discs — meteors, meteoroids — plasmas — turbulence 



1 INTRODUCTION 

Planets are formed in protoplanetary discs from interstellar 
dust. The electric charge state of the dust aggregates in the 
protoplanetary discs is one of the key parameters in under- 
standing a number of aspects of protoplanetary discs and 
protoplanetary formation. 

Planet formation begins with mutual sticking of /im- 
sized dust, most probably leadin g to extremely lo w density, 
fluffy structure of the dust (e. g. lOssenkopj [l993l ) . The oc- 
currence o f fluffy dust is sugge s ted by laborator y experi- 
ments (e. g IWurm fc Blumlll998l: iBlum et al.lll998h. by the- 
ories (e.g. lOrmel et al.ll2007l: IZsom fc pullemondll2008l). by 
N-body simulations (e.g. Kempf et al. 19991 : Suyama et al.l 
l2008l : IWada et al.l l2"008al ) . and by observations, including 
optical observatio ns of dust in star forming region (e.g. 
lEvans et al.ll200ll ) and observations of dust linear polariza- 
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tion in comet comae l|Levasseur-Regourd et al.ll2007h . For 



review of this field, see e.g. 



BlumI (|2"004D . 

The inner structure of the dust aggregates, relative ve- 
locity, and electric charge are key parameters that determine 
the growth a nd migration of du st aggregates. Dust rela- 
tive velocity l|Brauer et al. 20081) includes random motion 
caused by turbulence |Ormel fc Cuzzi|[2007h and Brownian 
motion (|Blum et al.lll996|). and bulk motion caus ed by verti- 
cal sedimentation IIDuUemond fc DominikllioMi ) and radial 
migration (IWeidenschilling||l977 ) . The col lision velocity gov- 
erns the growth rate ( Suyama et al.ll2008|'). compactificatio n 



IWeidling et al.ll2009l ). and disruption (|Wada et al.ll2008al ). 
of the dust. 

[Okuzuml Hooi) considered the charge state of the dust 
aggregates in protoplanetary discs. They assumed that the 
charge state is determined by absorption equilibrium of ions 
and free electrons. Since electrons have much larger thermal 
velocity compared to positive ions, plasma absorption makes 
all dust to charge weakly negative. The repulsive Coulomb 
force may suppress dust-dust collisional growth for all but 
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the heaviest dust species who can overcome the Coulomb 
barrier. 

It is also possible that the charge state of the dust is 
affected by dust-dust collision. The effect has been simply 
ignored in most research due to the fact that in protoplane- 
tary discs, dust has low number density, and is surrounded 
by weakly ionised plasma. We give quantitative estimates 
of the effect of dust-dust collision as a function of dust size, 
fractal dimension, and number density, and show for the first 
time that the dust-dust collision is actually an important 
factor in high dust number density regions of protoplane- 
tary discs. 

One of the possible dust-dust coUisional charging mech- 
anisms is known as the triboelectric process, where two bod- 
ies exchange electrons and sometimes molecular ions when 
they come into contact (e.g. lSickafoose et al.ll200ll ). Another 
mechanism is possible for materials with spo ntaneous sur- 
face charge, such as H2O ice crystals (e.g. iKudin fc Carl 
I2OO8I ). In this mechanis m the sur face matter within typical 



to electric discharge, they generally fragm ent without being 
thermally processed l|Guttler et ahlboosl ). 

Lightning in protoplanetary discs is strongly related 
to turbulence. The relative random velocity between the 
charged dust species that sets the dust to collide, results 
from the turbulence. Also the difference of the bulk velocities 
between the charged dust species that leads to macroscopic 
charge separation results from the turbulence. 

The turbulent state of the accretion discs is often ex- 
pressed in terms o f visco us a parameter introduced by 
IShakura fc SunvaevI (|l973l ). Since the specific angular mo- 
mentum increases outward in Keplerian discs, they satisfy 
Rayleigh's hydrodynamical stability criterion, and there are 
no clear mechan ism for hydrodyn amic turbulence in proto- 
planetary discs (|Sano et al.ll2004 ). On the other hand, the 
angular velocity decreases outward in Keplerian discs, they 
satisfy criterion for magnetorotational instability (MRI). 

Therefore, if a protoplanetary disc is ionised enough 
to sustain magnetic field, MHD turbulence is excited and 



depth ~ 1.0 X 10-* cm llMason fc Dash.200Q ) is exchanged Q parameter can be as large as 1.0 x 10 



1.0 X 10' 



together with contained charge. 

Surface space charge due to electron s pill-out is widel y 
known among metals and semiconductors (|Somoriailll994 ). 
the charge separation being ~ 10"^ cm deep for metals and 
~ 10"^ cm deep for semiconductors. H2O is unique in that 
molecular ions OH" and H3O''' holds the charge, and that 
proton e xchange betwe en the molecules (Grotthus mecha- 
nism; c.f. lAgmonl|l995l ) allows charge diffusion much faster 
than molecular ion diffusion. Thus su rface charge separ ation 
develops as deep as ~ 2.0 x 10"* cm (|Dash_e^^ 200j)^_Jbr 
example, ammonia lacks the mechanism i Goncalves et al.1 
1 19991 ). It is important that charge separation layer is deeper 
than exchange depth, because if the entire charge separation 
layer is exchanged, charge transport is neutral and coUisional 
charging do not take place. 

The dust-dust coUisional charging due to the exchange 
of this spontaneous surface charge of ice crystals, i s an estab- 
lished model in th e cont e xt of meteorology (e.g. iTakahasliil 
1 19781 : iBaker et al] Il987l : iDash et al.1 l20o!j | that explains 
lightning on earth. When two ice dust of different surface 
states collide, they exchange their surface charge, producing 
charged dust. When the charged particles within noncon- 
ducting gas are separated by some external force, electric 
field grows between them. At the point the electric field 
is larger than the dielectric field strength of the gas, rapid 
ionisation of the gas occurs, converting the electrostatic en- 
ergy into kinetic energy of the electrons and ions. This is 
electric discharge. Lightning in the earth's atmosphere is 
one of the most prominent, and well studied examples of 
electric discharge phenomena; in thunderclouds, typically 
3.0 X 10^" esu, or 1.0 x 10^ C of electric charge is repeat- 
edly separate d and neutralized with typical length scales 
1.0 x 10^ cm (jKoshak fc Krideij|l989l ). 

In protoplanetary discs, lightning is one of the can- 
didate mechanisms for chondrule heating, although com- 
pared to ot her models e.g. heating by shock wave (e.g. 
iMiura et al. I [2008). some difficulties have been pointed out 
i Weidnschillind Il997h . l^'or example, electric field cannot 
grow large enough to cause elec trostatic breakdown in stan- 
dard discs (|Gibbard et al.|[l997l ). Moreover, when mm-sized 
silicate aggregates made of /^m-sized monomers are subject 



ISano et al.l 1 199^ ). If the ionisation is suppressed, on the 
other hand, a ~ 1.0 x 10~^. For a typical protoplanetary 
disc it is believed that so-called 'dead zones' form between 
0.1 AU and 10 AU where inst abilities are da mped and 
gas flow is almost laminar (e.g. iGammij 1 19961 ) . But it is 
possible that MRI is active in the whole disc, if sufficient 
ionisatio n degree is mainta ined, for example by turbulent 
mixing llTurner et al.1 12007| ) or by self-sustained ionisation 
(|lnutsuka fc Sanol 20051 '). Thus ionisation state of the proto- 
planetary discs is critical in determining a and understand- 
ing the fate of planetesimals and protoplanetary discs (e.g. 
iKretke fc Linll2007l : iBrauer et al.|[200a l. 

The purpose of this paper is twofold: One is to solve 
the local charge exchange equilibrium of gas and dust nu- 
merically, for various dust parameters such as radii, fractal 
dimensions and dust number density, with dust-dust coUi- 
sional charging taken into consideration; Given the results, 
the other goal is to determine the critical dust number den- 
sity 77crit under which lightning to take place, as analytical 
functions of other dust parameters such as radii, fractal di- 
mensions and disc environment parameters such as temper- 
ature and gas number density. 

This paper is organized as follows. We define the terms 
we use in Table [T] and we list the symbols we frequently use 
in Table [21 In !}21we introduce the dynamic charge exchange 
equations and its equilibrium solution in schematic forms. 
We introduce circuit diagram to depict them (Fig. [l|. In SjS] 
we examine the processes in protoplanetary discs that set 
the parameters for the charge equilibrium equations. Crucial 
parameters are dust number density, the amount of charge 
exchange in single dust-dust collision, and relative velocity. 
In ^we estimate the electrostatic field strength, and define 
the critical number density r^crit for lightning in the proto- 
planetary discs. At this point all the equations are specified, 
and we solve them numerically. In ^ we show the results 
of the simulations. We describe four distinct phases of the 
charge distribution and explain the results using circuit di- 
agrams. We also give analytical estimates for electric field 
strength in protoplanetary discs and critical number density 
rj for lightning to occur. In JHJwe discuss the possibility of 
various phenomena caused by the highly charged dust and 
lightning in protoplanetary discs, and their observations. 
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particle 

gas dust 
plasma 

neutral gas ion electron smaller dust larger dust 

II I X I 

neutral cation anion cationic anionic 

Table 1. The terminology we use in this paper. 'Particle' is generic term for all components in the protoplanetary discs. Solid components 
are 'dust,' and the others are 'gas.' 'Gas' components are further subdivided into 'neutral gas,' and charged components, or 'plasma.' 
Finally, 'plasma' consists of 'electron,' the negative charge carrier, and various molecular 'ion,' the positive charge carrier. On the right 
side of the table, 'dust' is classified by their size as 'smaller dust' and 'larger dust.' Either can be 'anionic' or 'cationic' dust, depending 
on the material they consist of. We also use the one-letter symbols 'g', 'e', 'i', 'S', and 'L' for neutral gas, electron, ion. Smaller and 
Larger dust. The symbols for 'Cationic' and 'Anionic' dust are 'C and 'A'. We use variable I to represent one of these symbols. 



symbol 


value/dimension 


meaning 


definition 


r 


2.7 AU 


— constants — 

orbital radius considered 


- 




3.8 X 10^ g cm-3 


gas surface density of MMSN 






1.6 X 10-1 AU 


scale height of MMSN 






1.7 X 10^ K 


temperature of MMSN 


m 


1X1 


1.6 X 10-1" g cm-3 


gas density of MMSN 




' s 


1.6 X 10-12 g cm-3 


spatial density of smaller dust in MMSN 






1.6 X 10-" g cm-3 


spatial density of larger dust in MMSN 




Vch 


0.1 


charge exchange efficiency 


^3.4.11 


(Teh 


6.2 X lO'' e cm-2 


charge surface density 


^3.4.11 




3.4 X 10^ cm sec-i 


bulk velocity of larger dust to other species 








1 cLlHAUlli VciUCILV Oi pal LiClcB Ol BUCClco X 




All 


Q /I ■ — ^ 1 o3 r'ni car- 1 

o.^ X iu cm sec 


mean collision velocity between a smaller dust and a larger dust 




V 


1 


— independent variables — 

dust number density of the considered region 
divided by that of the MMSN model 




ri 


cm 


radius of a dust aggregate of species I 


p8)l 


Di 


1 


fractal dimension of a dust aggregate of species I 


|[28t 


mi 


cm 


— dependent variables — 

mass of a dust aggregate of species I 


l(29t 


Ps 


g cm— ^ 


condensed density of smaller dust 


VP^ 


Pl 


g cm— ^ 


condensed density of larger dust 




ni 


cm 


number density of dust of species I in condensed regions 


Pi /mi 


?i 


esu 


The charge carried by a single particle of species I 




Qi 


esu cm— 


The charge density carried by species I 


qiUi 


•^1,1' 


esu cm-^ s— 1 


charge transferred from species I to 
species I' per unit time per unit volume 


I144I49II 


'S'kiss 


9 

cm^ 


contact surface area within a dust-dust collision 


ll38t 




esu 


amount of charge exchanged within a dust-dust collision 


Vch ^ch '^kiss 


^COU 


9 

cm^ 


cross section between two charged particles 


ll30t.(l3Tl 


Jd 


esu cm-2 s— 1 


current carried by dust particles 




3p 


esu cm-2 s— 1 


current carried by plasma particles 


m 


Edis 


G 


critical electric field strength for lightning 


l56t 




G 


(local maximum of) electric field 





generated in the protoplanetary disc l|61|l 

X 1 whether the collision cross section between smaller dust and 

plasma particles are geometric(x «. 1) or Coulomb (x ^ 1) II113I I 



7?crit 1 the dust number density at which lightning takes place (1122 l l and II124I126I I 



Table 2. The list of symbols frequently used in this paper. 
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2 MODEL DESCRIPTION 

In this section we describe our models. In l2.1l we model the 
disc and the dust at the unperturbed state, then introduce 
the models for dust number density. In 12.21 we model the 
charge density and charge separation processes. 



2.1 Disc Model 

Unless otherwise mentioned, we focus on a local, uniform 
box at certain orbital radius r near the equatorial plane of 
the protoplanetary disc. We model the protoplanetary disc 
based on the m inimum-mass solar nebula (MMSN) model 
l|Havashi|[l98ll ). The gas surface density E^(r), disc scale 
height h!^{r) , and the temperature r'^(r) of the disc are 



E^(r) = 1.7 X 10- 



g cm 



ft-(r) =4.7x10-^ (^)' AU, 
T-(r) = 2.8xl0^(^)"* K, 



(1) 
(2) 
(3) 



where r is the distance from the central star. This leads to 
gas density distribution 



= 2.4xl0-«(-i^)"^ gcm- 



(4) 



The dust-to-gas ratio in M MSN is approxim a tely 1 .0 x 10 ^. 

We use the model by ICuzzi fc Zahnld l|2004 ). and m- 
troduce two species of dust, the smaller dust and the larger 
dust (see Table [1]) We further assume that surface density 
of the larger dust is 10 per cent of the total dust surface 
density. These two species are also either 'cationic' and 'an- 
ionic' The 'cationic' species receives the positive electric 
charge through dust-dust collision. See Appendix |A] for the 
justification of this two-dust model. We can also represent 
the role of vario us molecular ions by one abstract ion species 
'i,' according to lOkuzum il(|2009l'l. 

The motivation for this two-dust model is twofold. First, 
the two dust model is the simplest model that can handle 
the dust-dust coUisional charge separation and the macro- 
scopic relative velocity between the dust species. Second, 
the charge tendency of the dust and their size are strongly 
correlated. In one scenario, older dust are larger and also 
anionic. In another scenario, dust made of ice is larger and 
also cationic compared to dust made of silicate, (see H3.4l for 
the details.) Therefore, we expect that instead of consider- 
ing four (cationic smaller dust, cationic larger dust, anionic 
smaller dust, and anionic larger dust) species of dust, we can 
correlate the two size species with the two charge tendency 
species, (Table [T|, although both correspondences (smaller 
dust is cationic / larger dust is cationic) are possible. 

To summarise, we define the reference density of the 
smaller dust p'^ir) and the density of the larger dust p'^{r) 
as 



(r) = 1.0 X 10- V, (r) , 
(r) = 1.0 X 10-3p, (r) . 



(5) 
(6) 




Figure 1. The circuit diagram of the charge exchange process 
in dust plasma. Each arrow represents 'current' density J, which 
has the unit esucm~'^s~^, the amount of charge passed from one 
component to the other per unit disc volume per unit time. The 
arrow points from the component that receives negative charge to 
the component that receives positive charge. In this figure, i and 
e are ions and electrons created from ionising neutral disc gas. C 
and A are cationic and anionic dust defined in CAI 

Jei represents gas ionisation as 'current' from e vertex to i 
vertex; JiAi JiC^ JAe, and Jce are ion and electron absorption 

(n) 

to dust; Jac is dust-dust collisional charge separation and 
is neutralization current of charged dust-dust absorption. 

Alternatively, we can think of protoplanetary discs with dif- 
ferent gas or dust density than MMSN. We denote the ra- 
tio of the density of gas, smaller dust, and larger dust by 
fJgjVs 'Vl' respectively. Then the density of gas, smaller dust 
and larger dust is given by 



Pg = ngPg [r) , 
Ps = VsPs (^)' 

Pl = VlPl in ■ 

Mass of the smaller dust and the larger dust are and 
rrij^ , respectively. The number density is density divided by 
dust mass: 

n - 

m. 



(7) 
(8) 
(9) 



(10) 
(11) 



We further assume that within a local condensation re- 
gion, density for each component of the disc are multiplied. 



We estimate the mass as a function of the dust radius and 
the fractal dimension in ^3.11 



2.2 Charge exchange equations 

There are four species of charge carrier in our model — 
ions, electrons, cationic, and anionic dust (Table [l}. Charge 
exchange processes between these species are ionisation, 
plasma absorption, and dust-dust collision. The ionisation 
of the neutral gas molecules generates the ions and the elec- 
trons. Plasma absorption decreases the number of plasma 
particles and passes the lost charge to the dust aggregates. 
The dust aggregates also get charged by dust-dust collision. 

We label the particle species with letter I. The charge 
density carried by species I is Qi (the unit is esu cm"'^), and 
the charge transferred from species I to species I' is Ji ^/ (the 
unit is esu cm"'^ s"^). 

The charge density Qi of a species I is the product of 
their number density rii and their average charge per particle 
gi. For dust species, we assume that rii is known from num- 
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ber density model while qi is unknown; for ion and electrons 
we know qi but do not know rii. This constitutes the four dy- 
namical equations for four unknown variables q^,q^,ni, Ue 



^ = ^ Ja,C + J^,A - JA.e + fc,\) . (12) 

f^^i-(i.,c + Ac-ic.-jS) (13) 

^ = - {Js,^ - Ji,A - J^.C) (14) 

at qi 

^ = - (- Je,z + Ja,. + Jc,e) . (15) 
at (Je 

The current terms Ji ,/ are 

Ja,c = ^<lA,c'n'Anc'^A,c^v^,c (16) 

= \<ls\nAnc^^A^^^A,c (17) 

Ji^^ = enin_^(Tcon e) Ui (18) 

Ji,c = enin^CTcou (gc)e)i'i (19) 

J^,e = enen^cTcou (g^, -e)ue (20) 

Jp,e = enen^acouiqc,—e)ve (21) 

Je,« = Cng (22) 



where we have included neutral gas ionisation Je,i, dust- 
plasma absorption JA,i, Jc,i, JA,e, Jc,e, dust-dust colli- 
sional charge-up Ja,c, and dust-dust collisional neutraliza- 
tion terms. Here, Vi and Ve are the thermal velocity 
of the ions and the electrons, Ug is the number density of 
the neutral gas, is the ionisation rate, which is dominated 
by cosm ic ray ionisation near equat orial, r — 2.7 AU of 
MMSN (|Umebavashi fc Nakanol 120091 '). The exact value for 
these terms are given in 33] We have neglected, for example, 
the gas-phase recombination. 

We want to solve the equilibrium equations for the dy- 



namic equations (|12)I15[) : 

- Ja,c + Aa - JA.e + fc.\ = (23) 

Ja,c + Ac ~ Jc,e - fc.\ = (24) 

Je,i — Ji.A — Ji,C = (25) 

- Je,. + J A,. + Jc,. = 0, (26) 

together with charge neutrality equation: 

QA + Qc + Qr + Qe = 0. (27) 



We use circuit diagram (Fig. [T| to depict the dynami- 
cal equations p2ll5p . and to interpret the numerical equi- 
librium solutions (|23l27p in ^ The circuit diagram repre- 
sents charge-exchange processes; each vertex represents the 
species of charge reservoir and each arrow represents the 
charge exchange process. The size of the vertex circles rep- 
resents the amount of charge Qi. The thickness of the arrows 
represents the amount of charge transfer Ji,i'. We define the 
direction of the arrows so that the arrows point to the pos- 
itive charge receivers. 

In the system of equations depicted by a circuit dia- 
gram, charge density of each vertex Qi corresponds to an un- 
known quantities. Therefore, the number of unknown quan- 
tities is equal to the number of vertexes Ny- On the other 
hand, at the equilibrium, sum of the current flowing into 
each vertex is required to be zero (Kirchhoff's Laws); this 



gives us Nv equations but only A'^v' — 1 of them are indepen- 
dent. Charge neutrality gives us 1 equation. Thus we have 
A''v equations for Ny unknown values. 



3 CHARGE EQUILIBRIUM OF GAS AND 
DUST 

In this section we specify the current terms of the dynamic 
equations p6ll22p . especially the dust-dust collisional charg- 
ing terms Ja,c — ^^"a, by modelling the dust number density, 
structure, collisional cross section, surface charge exchange, 
and relative velocity. 

3.1 Fluffy dust model 

We use model of dust aggregates by IWada et al.1 (|2008bh . 
We consider dust aggregates composed of a large number 
of spherical monomers with radius =0.1 /^m. Each dust 
species I has its mass mi, the number of monomers that con- 
stitute the dust A'':, and representative radius ri. We define 
the fractal dimension of the fluffy dust Di in the following 
simple manner: 

7V:=(^)"\ (28) 

The dust mass is expressed in terms of monomer mass 
as follows: 

mi = mmNi — mm ( — 1 • (29) 

IWada et all (12008^ ) studies the coUision of the fluffy 
dust of the radii 1.0 x 10"^ ^ 9.1 x 10""* cm. The effect 
of offset collisions, collision between dust of much different 
sizes, and dust much larger than 9.1 x 10~* cm are yet to 
be confirmed. Therefore we make the following assumptions 
on smaller dust-larger dust collision. 

• If the smaller dust graze at the larger dust, i.e. if the 
line that passes the gravitational centre of the smaller dust 
and is parallel to the relative velocity vector do not intersect 
with the larger dust, the two dust aggregates do not stick to 
each other. Therefore the grazing cross section is of the order 
of r^Vj^. In this case they separate Ag^,^ of charge, which 
is the product of charge surface density ach and contact 
surface area Skiss • This contributes to the dust-dust charging 
current, Jc.a- 

• If the smaller dust bump into the larger dust, i.e. if the 
line that passes the gravitational centre of the smaller dust 
and is parallel to the relative velocity vector do intersect 
with the larger dust, the smaller dust do not penetrate the 
larger dust but becomes a part of the larger dust. The cross 
section is of the order Tj^^. In this case all the charges the 
smaller dust have are removed from the smaller dust charge 
density and added up to the larger dust charge density. This 
contributes to the dust-dust neutralization current, J^\- 

3.2 Collisional cross section of charged spherical 
object 

In this section, we estimate collisional cross sections for dust. 
The collisional cross sections for two electrically charged 



6 Takayuki Muranushi 




Figure 2. A grazing collision between a smaller dust (the blue 
solid sphere) and a larger dust (the red wire-frame sphere). The 
smaller dust creates a trench on the larger dust (the black cylin- 
der). 

spherical particle is given by 

(Tcou {q) = T^a exp {^-^^f^ {<ll' > O) , (30) 

(Tcou {q) = 7ra^ ^1 - ^J^^ {<ll' < O) (31) 

where q,q' is each particle's charge, T is the temperature of 
their relative mot ion and Tva is the geometric cross section 
fe.g. lSpitzej|l94ll ). 

Equation (I31|l represents the effect of Coulomb focus- 
ing: particles of the opposite charge attract each other and 
collide more often than when they are neutral. On the limit 
|gg'a~^| y=> ksT we can approximate the cross section as 
CTcou{q) — —TTaqq' (kBT)^^ , which is bi-linear on q and q' . 
On the other hand, cross section (|30|) represents the effect 
of Coulomb repulsion: for the collision between particles of 
the same charge only a portion of particles that belongs to 
the long tail of Boltzmann's distribution for temperature T 
can overcome the Coulomb barrier and collide. On the limit 
qq'a~^ ^ ksT the cross section vanishes quickly, but never 
reaches 0. 

We use Coulomb cross sections (|30|) . (|3H) to estimate 
the event rate of gas-dust collision and dust-dust collision. 

3.3 Collisional cross section and contact surface of 
fluffy dust 

The amount of charge exchanged in a collision, Ag^,^, is 
product of area of contact Skiss, upper limit of charge ex- 
changed per unit surface area of contact ach, and the non- 
dimensional efficiency factor rich- 

We leave the detailed argument to determine rjchCTch to 
>I3.4I Here we assume that r/chCch is known and describe 
how to estimate contact surface area Skiss- Since it requires 
another detailed simulation to estimate Skiss qualitatively. 



we resort to an order-of-magnitude estimate for this part of 
the work. 

We illustrate the collision between a smaller dust and 
a larger dust in Fig. [2] The smaller dust grazes the larger 
dust, pushes away the monomers that belong to the larger 
dust and creates a trench on the larger dust. The trench is 
a portion of the black cylinder in the figure. The radius and 
the length of the cylinder is and {r^r^^y^'^, respectively. 
Therefore, the surface area of the trench Sc is of order 

Sc^^//V,^/^ (32) 

and the number of monomers Nc required to fill the surface 
of the trench is 

iVc-r^-V/^V^^/^ (33) 

Their total surface area is also of the order of Sc- 

However, Skiss — ^s^^^'''l^^^ overestimates the actual 

contact surface area if the large dust is so fluffy that there 

is not enough monomers in the trenched volume to fill the 

trench surface. 

From the definition of the fractal dimension (I28|l . the 

number density of monomers within the larger dust material 

is 

(Af) ,r -3 -D, D, -3 /oA\ 

rv^ ' — Nr^ — rm ^r^^ ^ . (34) 
On the other hand the volume of the trench is 
W^r//V,^/^ (35) 
Therfore, the number of particle contained in the trench is 

f.r (A/)i/ ^ -D, 5/2 D. -5/2 ,na\ 

Nf = nl 'Vf Tm ^ r-g ' i ' , (36) 
and their total surface area is 

Q ^ 2 ,r ^ 2-D, 5/2 D, -5/2 

bF Tm Nf rm ^ r ^ ' r ^ ' . (37) 

If Nf < Nc, the surface of the trench is only partially 
covered by the monomers, and we estimate Skiss — Sf- On 
the other hand, if Nf > Nc, Np monomers are crushed 
onto the trenched surface, and since they overlap, about Nc 
monomers will take part in the charge exchange. In this case 
we estimate Skiss — Sc- To summarize, we assume that Skiss 
is the smaller of ((32)1 or (p7)) : 

o • f 3/2 1/2 2-13, 5/2 -5/2\ /oo^ 

Skiss = mm (r^ ' ' , '"s '^r, ^ j • (38) 

3.4 Charge separation processes 

There are generally two classes of possible charge separation 
processes in protoplanetary discs. 

One is surface charge exchange, where each dust has 
some kind of spontaneous charge separation (|Kudin fc Carl 
|2008| ). so at the initial condition each dust charge is zero as 
a whole (globally neutral), but there are charge separation 
within the dust particles (locally charged) . For example, wa- 
ter ice crystals tend to gather negative charge at its surface 
and positive charge inside. When two dust aggregates with 
different charge collide and melt partially, they exchange 
molten material and the charge included in the molten ma- 
terial. As a result each dust gets globally charged. 

The other charge separation mechanis m may be tribo- 
electric processes (e.g. iDesch fc Cu^l2000D . In this 
the initial condition each dust is both globally and locally 
neutral. When two dust aggregates made of materials with 
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different electron affinity collide, the surface electrons move 
from one material to the other. As a result each dust gets 
globally charged. 



3.4-1 Surface charge exchange I — larger dust is anionic 

The mechanism we consider the most plausible for the dust- 
dust collisional charge separation is surface charge exchange 
betwe en ice dust. For t he du st aggregate of ice mantled sil- 
icate, ICuzzi fc Zahnlel (|2004 ) proposed a condensation sce- 
nario, that at the snow line ice larger dust drifting inward 
dissociate and many smaller dust form. 

There are established models on charge separation 
caused by ice-ice dust collision in the context of t hunde r- 
cloud meteorology (for review, see e.g. iDash et al.l l|200j) '). 
We will carefully import them as a charge separation model 
in protoplanetary discs. The essential steps to cause light- 
ning on earth are (1) spontaneous charge separations on ice 
crystal surfaces, (2) existence of different dust species with 
different spontaneous charge separation per surface area, (3) 
collisions between the different dust that leads to global 
charging of each dust and (4) relative motion between the 
globally charged dust to create electrostatic field. 

For (1), we argue that the charge separation per sur- 
face area is quantitatively the same as the values measured 
in laboratory experiments. For (2), dominating dust species 
in charge separation process in protoplanetary discs is un- 
certain, and we discuss two possibilities (c.f. i]3.4.1l ^3.4.21 ) 
in this work. For (3) and (4), we make simple estimations 
for the collision rate and relative velocity in protoplanetary 
discs. 

Ice crystal surface is intrinsically charge-separated. Ice 
is negatively charged near the surface, and the inside is posi- 
tive. The typical charge surface density for stable ice surface 
is (Jcfi — 3.0 esu cm"'^ or a^h — 6.2 x 10® e cm~^ and the typ- 
ical skin depth of the charged layer is dch — 2.0 x 10~^ cm, 
though charge surface d ensity for fast-gr owing ice surfaces 
are larger and shallower (|Dash et al.ll200i] ) . This charge sep- 
aration has a general explanation as a result of interaction 
between hydroxide(0_g ~) and hydronium ( H^O^) ions and 
a hydrophobic surface l|Kudin fc Caij [20081 ') . and the above 
value of typical charge surface density is observed at liq- 
uid w ater-air surfaces as well as at ice crystal-air surfaces 
((Takahashi 2005 ). Therefore we use the value for ice- vacuum 
surfaces as well. 

In the thundercloud, there are varieties of ice crystals 
with different surface charge densities, depending on the sur- 
face history of the ice crystals. Newly formed surfaces have 
larger charge surface density than old surfaces, because they 
have higher fractal dimension and deeper amorphous layers. 

We now co nsider how surface ch arge exchange works 
in the model of ICuzzi fc Zahnlel ((2004 1. Larger dust that 
migrate towards the snow line has old surface and has less 
negative charge surface density, while smaller dust formed 
at the snow line have new surface and larger negative charge 
surface density, as in meteorological case. Note that before 
collision each dust is globally neutral. 

At the collision, the surface of the dust aggregates melts 
and the surface charge density is exchanged, and averaged. 
The larger dust, having less surface charge density than the 
smaller dust, receives more negative charge than it gives. 



Therefore the larger dust becomes anionic, smaller dust be- 
comes cationic. 

Laboratory experiments (lTa.kahashi"l978'), in-situ ob- 
servations and meteo rological estimates (,Gaskell et al.lll97§ : 
IChristian et al.lllQSOl ') suggest that for mm-size ice crystals, 
at least 10 per cent of the total surface charge within con- 
tact surface is e x chang e d in a single c ollisio n; experiments by 
iMason fc DashI (|2000l ): iDash et all (120011 ) suggests almost 
rich = 1.0. As a conservative estimate, we use rich = 0.1 
unless mentioned otherwise. 



3.4-2 Surface charge exchange II — larger dust is cationic 

It may be possible that charge separation processes occur- 
ring in protoplanetary discs are different from those occur- 
ring in the terrestrial thunderclouds. The collision time-scale 
in the protoplanetary discs is much longer than th at in a 
thun dercloud, so long that sintering may take place (|Sirond 
Il999iy As a result. The surface state of old ice larger dust 
and young ice smaller dust might resemble each other. If 
they are identical, some random charge exchange by colli- 
sion is still possible, but they do not exchange charge on 
average. 

However, compared to thundercloud, protoplanetary 
discs are more dirty and fine-grained; they contain much 
dust made of materials other than ice such as silicates, and 
the monomer size is 0.1 /im rather than 1 mm. Since the 
monomer size is smaller than typical skin depth of the charge 
separation dch — 2.0 x 10~* cm mentioned above, it is possi- 
ble that ice smaller dust and silicate smaller dust with thin 
ice mantles formed at the snow line is inefficient in separat- 
ing charge. There may be silicate aggregates with no sur- 
face charge separation. Meanwhile old larger dust that have 
travelled from the far end of the protoplanetary disc have 
undergone sintering and have developed thick mantles with 
full surface charge separation. 

In such scenario, the larger dust has more surface charge 
separation than the smaller dust. Therefore, collision be- 
tween a larger dust and a smaller dust still leads to charge 
separation but the larger dust becomes cationic, and the 
smaller dust is anionic in this case. We assume that r/ch = 0.1 
and (Jcfi — —3.0 esu cm~^ in this case (The charge exchange 
rate has the same magnitude but the opposite sign compared 
to that of EXU) 

Both scenarios, the larger dust is anionic and the larger 
dust is cationic are plausible. They may even take place in 
the different parts of the same disc simultaneously. There- 
fore, we have decided to take both scenarios into consid- 
eration. To that end, we treat the concept of cationic and 
anionic dust separately from the size of the dust. 



3-4-3 Trihoelectric charge separation 

iDesch fc Cu^ (|2000'1 have proposed that collision between 
large silicate grains and fine iron metal grains leads to tri- 
hoelectric charge separation. For instance, silicate dust of 
radius 3.0 x 10~^ cm will gain 5.4 x 10^ e charges per dust. 
The process can be built into our model in the same manner 
as we treat surface charge exchange processes. 
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D, m, St 



,,turb 



1.0 X 10^ cm 3.0 3.9 x lO** g 1.6 x 10^ 7.8 x 10^ cm s^^ 1.1 x 10^ cm s^^ 
1.0 cm 3.0 3.9 g 1.6 x 10"^ 2.1 x 10^ cm 8"^ 5.9 x 10^ cm s"^ 

1.0 x 10^ cm 2.4 1.5 X 10^ g 6.2 x IQ-"* 7.9 cm s"! 3.6 x 10^ cm s"! 

Table 3. T he estimated Stokes number, the bulk velocity due to the inward migration u™^ jBrauer et al.|[2008l) . and the turbulent 
speed u l|Ormel fc Cu2zill2007l 'l for some typical large dust parameters, for MMSN equatorial at r = 2.7 AU. 



3.5 Relative velocity 

When a cloud of positively and negatively charged dust is 
separated much larger than plasma Debye length 



Xd = 



5.0 X 10- 



(39) 



the electrostatic field between them become observable. In 
order to cause such macroscopic charge separation, there 
must be a significant relative bulk motion between anionic 
and cationic dust. Inward migration of large dust is a source 
of this bu lk motion. The sedime ntation may act in the same 
way. Also lDesch fc Cuzzil lj2000h have proposed that largest 
eddies in turbulence of protoplanetary discs cause bulk mo- 
tion between smaller dust and larger dust. Such effects on 
the relative velocity betwe en dust species in MMSN has been 
studied (see lBrauer et af] (|2008 ) and references therein). 

Here, we simply assume that the largest contribution 
to the smaller dust-larger dust relative velocity is the bulk 
motion of the larger dust, and the velocity is Au^^^g = = 
3.4 X 10^ cm s~^, the catastrophic collision ve locity of the 
ice du st aggregates of 9.1 x 10~* cm size dust l|Wada et al.l 
l2008bl ). Note that the non-sticking velocity threshold de - 
crease as the monomer size increase (|Blum fc Wurmll2000l ). 
We also check our analytic formulae with smaller values of 
Avj^^g and Uj^ assumed. 

Dust migration speed are comparable to this value at 
some stages of the dust growth. On the other hand, tur- 
bulent motion is faster than the value for most of our pa- 
rameter range (c.f. Table [S}. Turbulent mode that is larger 
than the scale of interest can be treated as bulk motion, and 
can be used to explain the charge separations of the scale. 
The scale can be as la rge as of order of disc scaleheight 
jBalbus fc Hawlevlll99lh . 



3.6 The charge equilibrium equations 

By substituting the results of analyses up to here into (|12I 
lisp we have the following dynamic equation for charge trans- 
port: 



dt 
dQg 

dt 
dQ, 

dt 
dQe 

dt 



(40) 
(41) 
(42) 
(43) 



where the current density terms p6ll22p become: 
/ 2rg . g„ \ 



, Ql Qs 1 

1 If r I L ^ 



^ cou 
~Q ^ cou 



Ql 



,e,r„,kBT Vi 



(44) 

(45) 

(46) 

(47) 

(48) 
(49) 



In (|44l) . the amount of current exchange Ag^^^ is prod- 
uct of contact surface area Skiss and surface charge density 
(Teh , each described in il3.3l and ^3.41 The contact surface 
area Skiss is the function of dust radii and dust fractal di- 
mensions; see equation (I38|l . The surface charge density cTch 
depends on the dust material. The relative velocity of the 
larger dust and the smaller dust is Au^^.g — 3.4 x 10"^ cms~^, 
as we have discussed in i]3.5l The cross section term (Jcou is 
the Coulomb cross section introduced in ti3.2l We assume 
Vi and Ve to be thermal velocities of ions and electrons. For 
ionisation in MMSN at r = 2.7 AU, cosmic ray ionisation is 
the m ain contributor and C, ~ lO"'^^ l|Umebavashi fc Nakanol 
l2009h . We introduce the nondimensional dust number den- 
sity r] (dust number density in unit of MMSN values), so 
that in equations (OU]), rig = 1, and Vs ~ Vl ~ V- From 
those density term, the number density terms ng,ng, are 
given as pg /rUg, pg /m^ , /m^. The masses of dust aggre- 



sions; see equation 

All the variables that appear in the current density 
terms (|44)I49[) are controlled by five parameters; radii of the 
dust aggregates (rg,r^), their fractal dimension {D^,D^), 
and the nondimensional dust number density r]. 

The equilibrium equations (I23II27|I become: 



Jr 



Jl 
J. 



J^,e = 
J..e = 




Ji.e 



(50) 
(51) 
(52) 
(53) 
(54) 



Again note that, out of four Kirchhoff's Laws (|50l53p only 
three of them are independent, and the charge neutrality 
condition (|54p is necessary. 
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4 CRITICAL DUST NUMBER DENSITY FOR 
LIGHTNING 

In this section we derive the strength of electric field gener- 
ated by the relative motion of the large and small dust, and 
set conditions for macroscopic electric discharge events, or 
lightning. 

Lightning occurs when the maximum electric field in 
the plasma -Emax exceeds the critical value -Edis- The crit- 
ical electric field Edis is determined by the condition that 
an electron accelerated by the field has kinetic energy large 
enough to ionise a neutral gas molecule. Let Zmfp be the 
mean free path for electron. Then an electron accelerated 
in electric field of strength E receive the energy of order 
eElratp. The ionisation potentials AWion for H, H2, and He 
molecules are 13.6 eV, 15.4 eV, and 24.6 eV respectively 
l|Dulev fc Williams! 1 19841 ). We use AW^ion = 15.4 eV in this 
work. Therefore the critical value _Edis of electric field for the 
lightning satisfies: 



AWion 

-fidis 



(55) 
(56) 



Next we derive the value of -Emax- When the differential 
motion between the oppositely charged dust species con- 
tinues much longer than the plasma Debye length, it can 
be interpreted as current carried by the dust jo generat- 
ing electrostatic field, and the plasma counter-current jp is 
induced in the neutralizing direction . We consider that jp 
is carried by electrons, and neglect current carried by pos- 
itive ions because it is at most the same order as that by 
electrons. Moreover, even if positive ions are accelerated to 
AlVion and ionise other molecules, they increase the electron 
number density only linearly, not exponentially. 

The dust current jo is estimated simply, by the product 
of dust charge density and macroscopic motion Uj^ , as: 



(57) 



On the other hand the particle current jp is determined 
by the Ohm's law: 



jp — ^^-Emax, 

where f is the electric conductivity. 



(58) 



(59) 



me Ve 

£max is determined at the equilibrium of these two currents 
jn and >: 

jD+jp = 0. (60) 
By substituting ((53, dlHJ, and ((59l) into ([eO]), we obtain 



Smax = - ""7^^:^ . (61) 

Now that we know both Bmax and E^is, the condition 
for electric discharge is 



|E^max| ^ E^dis- (62) 

By substituting (|56p and (|6H) into (|62p , we have the follow- 
ing form of the condition for electric discharge: 



Within our parameter range of interest, the behaviour 
of the left hand side of lf63)l as we increase r) is that it first 
keeps values much smaller than the right hand side and then 
it monotonically increases (c.f. Figure. [3] [4]). Thus there is 
a unique value of rj at which the equality for (|63|l holds. 
We define this value to be 77crit, the critical dust number 
density at which lightning takes place. Note that the condi- 
tion doesn't depend on the detail of the electron stopping 
processes because we can eliminate /mfp from the condition. 



5 RESULTS 

We have performed two sets of numerical experiments. In 
the first set of experiments, we fixed the set of parameters, 
Vg, r^, Dg , and to some typical values. We varied the 
dust number density r), and calculated charge density for 
each species of particles at the equilibrium. 

In the second set of numerical experiments, we varied 
the set of input parameters, r^, r^, D^, and D^, and for 
each set of input parameters we calculated the dust number 
density required to cause electric discharge r^crit. 

For all these simulations we assumed the environment at 
the equatorial plane and the snowline of the MMSN model; 
r = 2.7 AU, = 1.7 x 10^ K, = 1.6 x 10" 
= 1.6 X 10~" g cm-^ f. 

The results of the first set of experiments are in i]5.1l We 
found that the dust-plasma system experience four phases 
as we increase rj. We interpret this result in tj5.2l The results 
of the second set of experiments are in i)5.3l We derive the 
analytic formula for rjcrix in i]5.4l 



-1" g cm-^ 
1.6 X 10"^^ g cm"^. 



AW^io 



(63) 



5.1 Equilibrium charge density of particles as a 
function of dust number density 

We found that as we increase 77 while keeping other dust 
parameters constant, the equilibrium charge densities Qi — 
qiUj experience four phases (Table|4]). Fig.[3]and Fig.|4]shows 
the typical four phases behaviour. 

In this and the next sections, we explain the origin of 
the four phases, using the circuit diagrams (Fig.[S| as a great 
help. The four-phase behaviour we describe here is indepen- 
dent of most of the details of charge exchange processes. In 
fact Fig.Omodel and Fig.Umodel have the opposite sign for 
dust-dust collisional charge exchange, but the evolutions are 
almost similar. The rest of the discussion in following sec- 
tions is based on the former case, which we consider is most 
plausible (see ijSAl}. The discussion is easily generalized to 
the other case. 

To analyse the result, we first identify the dominant 
processes by comparing the competitive current in circuit 
diagram, then write down all the unknown values in simple 
polynomials of rj. Fig. [S] illustrates the transition of dom- 
inant process in the circuit as dust number density rj in- 
creases. The two particles with the largest charge density is 
marked by larger circle. There are always two of them, one 
carrying most of the system's positive charge and the other 
negative, thus charge neutrality holds. The arrows and their 
line width represents direction and amount of currents. La- 
bels for dominant currents are marked with thick rectangle, 
sub-dominant currents with thin rectangle, negligible cur- 
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(a) 


ion-electron 
plasma phase 


(in this paper 2 \ Qe\ > Qi) 


(b) 


ion-dust 
plasma phase 




2\Qe 

Ja,c 


\<Qi , 

1 < \JA,e\ 


(c) 


charge- up 
phase 




1 < •^A.C 1 ^ 1 •^'■tC 1 


(d) 


dust phase 






<\Ja,c\ 



Table 4. The names and conditions for four phases of charge sep- 
aration. They are basically named after dominant charge carrier 
of each phase. 



rents with dashed rectangle. The names and conditions for 
each phase is listed in Table |4] 

There are two major consequences of the size difference. 
Larger dust is much fewer in number density. So in the fewer 
dust limit (t; <<< 1) the larger dust carries much less charge 
density than smaller dust do. Since larger dust is the fewer, 
one larger dust collides with smaller dust much more often 
than one smaller dust does with larger dust. Therefore larger 
dust are the species that experience the quick charge density 
raise in (c)charge-up phase. The main role of the smaller 
dust is to absorb plasma and keep the charge neutrality. 

5.2 Four phases of charge separation as a function 
of dust number density 

5.2.1 Ion-electron plasma phase 

In ion-electron plasma phase (Fig. [5] (a)), the dominant path 
of charge transfer is 



e" ^ i ^ C ^ e~, 
the next-dominant path is 
-* e" . 

Therefore, we have following current hierarchy: 



» Ji, 
» 



-J. 



,0 



(64) 
(65) 

(66) 



The amount of current for path H64|) is constrained by edge 
e~ — > i^; since we have assumed that C, and rig is indepen- 
dent of Tj, so is Je,i. 

From charge neutrality (|27p . Qe = Qi and therefore 
Ue = Hi. So equation Ji,^ ~ J^^e is satisfied by setting, in 
equations (fT^ and (^1) . 



i^cou {qc ,e)vi = acou {Qc , — e) Ve, 
(Tcou (Qc 1 e) 



Ve 

— oc r; . 

Vi 



(67) 
(68) 



(Jcou {qc,-e) 

Equation tells us that crcou{qc , e) /ccouiQc , '^^i'^" 
stant of ri. This means oc rj^ because the only ?7-dependent 
term in acou is q^. . By definition of dust number density fac- 
tor 77, rip <x rj^ , so oc rj^ . 

By similar argument we can deduce oc rj^ from 

■^'■A — '^A,e- 

In other hand, to satisfy Ji,^ oc rj^ and ~ J^^e oc ?7° we 
need ni,ne oc ry"^. And since qi,qe oc r;", we have Qi,Qe oc 



In this phase, ions and electrons are the major carriers 
of positive and negative charge. Equation (|68)) also tells us 
that acou{qc,e)/acou{qc,-e) — Ve/vi ^ 1. This is inter- 
preted as follows: Since thermal velocity of electron is much 
faster than that of molecular ions, electron is more rapidly 
absorbed to neutral dust than ions. Therefore dust contin- 
ues to acquire negative charge, until its negative charge is 
enough to repulse most of the electrons inflow to attain 
a current equilibrium. Both cationic and anionic dust are 
forced to charge negative to hold back the overwhelming 
electron absorption. 

To summarise. 



Qi <xri 
Qe oc rj' 

Qa « n'^ 
Qc v 



+1 



(69) 
(70) 
(71) 
(72) 



5.2.2 Ion-dust plasma phase 



The system enters ion-dust plasma phase when the negative 
charge in dust Q^ become comparable to that in plasma Qe. 
Charge neutrality (|27[) requires free electrons to decrease. 
So the Coulomb barrier of dust species become weaker until 
Coulomb cross section approximates geometric cross section 
acouiqc)^) — '^cou{qc)^) — ""t'c^ °^ where electrons and 
ions are equally absorbed to the dust. 

In ion-dust plasma phase (Fig. [5] (b)), the dominant 
path is still 



^ i C — > e~, 
and the next-dominant path is still 

i+ — > A ^ e", 

and the same current hierarchy holds: 



Je.i - 
» J^,A 
» J^.n 



- J A .e 



Jr: 



oc rj 



(73) 



(74) 



(75) 



However, now that crcoti(?c,e) ~ cjcouiqc^^), equation 



Jp,e is satisfied by setting, in equations (flQll and 



UiVi = UeVe, (76) 

— ~ — OC r; . (77) 

rie Vi 

So the ratio Ui/ne is kept constant to Ve/vi = 6.1 x 10^. 
Still, in order to have Ji^^ oc rf and ~ J^.e oc nf we need 
Hi, Ue OC ri~^ . Since g^, qe oc r;", we have Qi, Qe oc j;^^. 

In this phase the cationic dust carry most of the neg- 
ative charge while ions carry most of the positive charge of 
the system. Therefore, the charge neutrality equation (|27|l is 
dominated by these two components, and Q^ oc Qi oc r]~^ . 

In this phase anionic dust also feels the same environ- 
ment as cationic dust, so oc rj~^ . However as iq ap- 
proaches to (c)charge-up phase, dust-dust coUisional charge 
separation J^,^ gradually comes into play and increases. 
Therefore in Fig. Owe can see the power law oc r)~^ only 
at the beginning of (b)ion-dust plasma phase. 

To summarise. 



Qi OCT] 



(78) 
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symbol 



value 



r 


2.7 AU 


vixi 


o.o X iU g cm 




1.6 X 10-1 AU 




1.7 X 10^ K 




1.6 X 10^1° g cm-^ 




1.6 X 10^12 g cm-3 


pr 


1.6 X 10-1^ g cm-3 




9.3 X 10-1 g cm-3 




1.0 X lO-"' cm 




1.0 X 10-* cm 


^I, 


1.0 X 10^ cm 




3.0 


Dr. 


3.0 


c 


1.0 X 10-1* 




3.4 X 103 cm s-i 




3.4 X 10^ cm s-i 




6.2 X 10^ e cm-2 


Vch 


1.0 X IQ-l 



Figure 3. Amount of charge stored in each species, erie, en^, \qg \ , and \ as functions of rj. This figure is for ice dust-ice dust 
case, so larger dust is anionic and smaller dust is cationic. The polarity matches that of Fig. |5] The radius of smaller dust, radius of 
larger dust, fractal dimension of smaller dust, fractal dimension of larger dust are 1.0 X IQ-* cm, 1.0 X lO'^ cm, 3.0, and 3.0 respectively, 
(a), (b), (c), and (d) corresponds to the four phases described in i|5.2l The yellow arrow denotes the critical number density rj where the 
macroscopic electric discharge condition II63I I is met. The settings of the simulation that produces this figure is in the right table. 




symbol 


value 


r 


2.7 AU 




3.8 X 10^ g cm-3 




1.6 X 10-1 AU 




1.7 X 102 K 




1.6 X IQ-io g cm-3 




1.6 X 10-12 g cm-3 




1.6 X 10-13 g cm-3 


Pm 


9.3 X 10-1 g cm-3 




1.0 X IQ-"^ cm 




1.0 X 10-1 cm 


rL 


1.0 X 10^ cm 


Ds 


3.0 


Dr. 


3.0 


c 


1.0 X 10-1* 




3.4 X 10^ cm s-i 




3.4 X 10^ cm s-l 


<^ch 


-6.2 X 10^ e cm-2 * 


Vch 


1.0 X 10-1 



Figure 4. Amount of charge stored in each species, erie, erti, \q^\ng, and |g^|n^, as functions of r]. This figure is for ice dust-silicate 
dust case, so the larger dust is cationic. Radii and fractal dimensions of dust, and other parameters are all same as in Fig.|3] except that 
the amount of charge exchanged in a collision has the opposite sign, so larger dust is cationic and smaller dust is anionic. 
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Qa « n~^, 



oc rj 



(79) 
(80) 
(81) 



In ion-electron plasma phase and ion-dust plasma phase 
the dust-dust coUisional charging is ineffective. So we can 
understand these two phase without dust-dust coUisional 
charging (see iOkuzumil (20091 ) and references therein.) 



5.2.3 Charge-up phase 

The system enters (c)charge-up phase when J^.c becomes 
larger than J^,e- Now anionic dust has their own negative 
charge supply from dust-dust collision, their negative charge 
grow quickly, and acou{qA i ~e) become rapidly small. At this 
point, the circuit switches one of its current path. 

In charge- up phase (Fig. [5] (c)), the dominant path is 

stiU 



e ^ i ^ C ^ e , 

but the next-dominant path is 



■+ 



(82) 
(83) 



The amount of current for path H82|) is constrained by edge 
e~ — + i^; since we have assumed that ^ and rig is indepen- 
dent of rj, so is Je,i. The amount of current for path (|83p 
is constrained by edge A — > C (I16|) : since we have assumed 
that Ag^_p is independent of ri, Ja,c ■ 
Therefore, we have following hierarchy: 



» Ji,A ' 
» JA,e- 



Jr: 



(X rj 
oc ri^ 



(84) 



The path (|82|) is as same in ion-dust plasma phase, 
leading to Qi,Qe oc rj~^ , and charge neutrahty requires 
Qc ^ V~^- 

In dust charge-up phase, however, anionic dust has so 
much charge that electrostatic potential for electron and ion 
at the surface of larger dust is larger than their thermal 
energy; this is qq'a~^ >>> ksT limit of the Coulomb cross 
section ((30)) . (f30)) . Thus Gcou{qA , — e) and Gcou{qA , e) oc 
9a i'^ (|18|) . Substituting rii oc rj~^ and oc r}^ into Ji,^ oc 
rf , we hav 

To summarise. 



oc rf and « rf . 



(85) 
(86) 
(87) 



Qi oc ri~ 
Qe oc ri~ 

Qa V 
Qc <x rj 

At this phase, by substituting equations (|86p (|87p into 
equation (|61l) we have 



En 



oc rj 



(89) 



The -Bmax has the dependency of »7* in this phase, instead 
of i? o c Tj^ dependence used, for example, in lGibbard et all 
([l993). Moreover, at the end of dust charge-up phase there 
is a steep increase in and steep decrease in Qe. These 
means that the electric discharge condition (|63p meets at 
smaller value of rj. 



5.2.4 Dust phase 

The system enters (d)dust phase when J^,^ becomes larger 
than Ji.p . Now the charge states of both anionic and cationic 
dust is governed by dust-dust collision, and the plasma com- 
ponent is sub-dominant to the dust. 

In dust phase (Fig. [5] (d)), the dominant path is 



(90) 



the dust-dust collision is now short-circuiting. The next- 
dominant path is 



C 



A. 



(91) 



The amount of current for path (|90|) is constrained by edge 
A —* C Uni; since we have assumed that Ag^^^ is indepen- 
dent of 77, J^,c oc 77^. 

The amount of current for path (|91|) is constrained by 
edge e~ ^ i"^; since we have assumed that C, and rig is 
independent of 77, so is J^.i. 

Therefore, we have following hierarchy: 



J A>C 

» Jc,^ 

» JA,e- 



•^CA 



OC »7 
oc»7° 



(92) 



Q. 



Equation J^.^ ~ Jc,a requires Aq^^^'^A.c ~ 

Therefore only rj dependent term q^, must satisfy 
oc rj^' , leading to oc Tj^ . Charge neutrality leads to 
oc rj^ . 

The path (|9ip gives us Qi,Qe oc rj"^ , same as in ion- 
dust plasma phase and in dust charge-up phase. 

At the boundary of (c)charge-up phase and (d)dust 
phase there is a jump of dust charge. This is because when rj 
cross the boundary dust charge grows until dust-dust coUi- 
sional neutralization can compensate dust-dust charge sep- 
aration. 



To summarise, 



Q, 


oc 


-1 

V 


Qe 


oc 


-1 

V 








Qa 


oc 


V^' 








Qc 


oc 





(93) 
(94) 
(95) 
(96) 



5.3 Critical dust number density as function of 
dust parameters 

We now explain the details of the second numerical exper- 
iments, where we varied the set of input parameters, r^, 
, Dg , and , and for each set of input parameters we 
calculated the critical dust number density r/crit at which 
the lightning strikes. The numerical results strongly suggest 
that the parameter space (r^ , r^^ , , ) is subdivided into 
several regions, at each of which r^crit is a simple analytic 
function of parameters (r^ , , , D^). 
The parameter ranges are 



1.0 X 10"* cm < < 1.0 x 10^ cm, 

1.0 cm < < 1.0 X 10^ cm, 

2.0 < Ds < 3.0, 

2.0 < < 3.0, 



(97) 
(98) 
(99) 
(100) 
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with additional constraints 



(a) ion-electron plasma phase 




(c) charge-up phase 




Ds < D^, 
1.0 < 77 < 1.0 X 10'^ 



(d) dust phase 



(101) 
(102) 
(103) 

Constraint (llOlll requires that the smaller dust is smaller 
than the larger dust. Constraint (|102[) comes from em- 
pirical fact that larger dust aggregates have experienced 
more compactification, and have higher fractal dimension 
ISuvama et allbOOSl : IWada et al.li2008a l. Constraint ([10311 
is cutoff value of our computation. 

We visualize the four-dimensional field 
Vc-ritir g , r ^ , D g , D ^) in the figures at the last of this 
paper, by choosing some representative points and present- 
ing several 2-dimensional sections that passes the point. As 
the mass, the radius, and the fractal dimension of a dust is 
related by equation (I29|) . we have some freedom of choosing 
the direction of 2-dimensional section. We keep mi constant 
when we vary Di (the dust puff up with constant mass); we 
keep Di constant when we vary ri (the dust mass increase 
with constant fractal dimension). 

First, Fig.|6]shows the 'fluffy dust' cross sections, where 
the representative dust are = 1.0 x 10~^ cm, r^^ — 1.0 x 
10^ cm, = 3.9 x 10"^ g, = 1.5 x 10^ g, = 2.0, and 
Dj^ — 2.4. The critical number density is r^crit ~ 7.37 x 10^ 
for this representative parameter. 

The second set of Fig. [7] uses the 'hard dust' cross sec- 
tions, where — 1.0 x 10"* cm, = 1.0 cm, = 
1.9 X 10"^^ g, = 3.9 g, = 2.7, and = 3.0. The 
critical number density is r^crit = 3.01 X 10^ for this repre- 
sentative parameter. 

The third set of Fig. [S] is the r^crit averaged over the 
parameters that do not appear in the axes, to show the 
tendency of overall dependence on the parameters, and to 
demonstrate the precision of the analytic formulae. 

The fourth set of Fig. [9] uses the same representative 
dust as in Fig. [S] but is the result of another simulations, 
where we are now extremely pessimistic and assume that 
the charge exchange is four orders of magnitude inefficient 
iVch = 1.0 X 10~^ instead of rjch = 1.0 x 10"^). Even though, 
the number density j^crit required for lightning has raised 
only by two order of one magnitude. The critical number 
density is rjcrit = 6.59 x 10'' for the representative parameter. 

The fifth set of Fig. [TD] shows the averaged ?7crit for the 
pessimistic case rjch ~ 1.0 x 10~^. We later examine the 
accuracy of our formulae by fitting Fig. [TU] with the formulae 
using correction factors determined by Fig. [S] data. 



Figure 5. The evolution of the charge density and current density 
as dust condense. As dust number density rj increase, J^^i oc r?" 
is constant while J^,^ °^ v'^ grows, and the particle experience 
four phases in order (a) — > (b) — > (c) (d). (a) At ion-electron 
plasma phase, most of the charge is carried by plasma species and 
the charge state of the dust is governed by plasma absorption, (b) 
At ion-dust plasma phase, the current balances are same as it was 
in ion-electron plasma phase, but now the negative charge carrier 
is cationic smaller dust, (c) At charge-up phase, anionic larger 
dust has sufRcient charge to cut off J^,e. (d) At dust phase, most 
of the charge is carried by dust species and the collisional charging 



5.4 Analytic formulae for lightning conditions 

In this section we derive the analytic form, of T^crit and light- 
ning conditions. Numerical results obtained in i]5.3l are of 
great help in deriving these analytic formulae. We show at 
the end of i]5.4.4l that by our analytic formulae we can fit 
364325 numerically-obtained points distributed among six 
decades with 21 per cent precision. Moreover, the formulae 
'predicts' results of another simulation with 59 per cent pre- 
cision, where charge exchange is 10* times inefficient. These 
agreements are good evidences for correctness of both nu- 
merical and analytical results. 

We made plots like Fig. [3] and Fig. [4] for many points 
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within our parameter space, and found that r) = Tjcrit is met 
at the boundary of (c)charge-up phase and (d)dust phase in 
most cases, and sometimes in (d)dust phase or (c)charge-up 
phase. Therefore we derive analytic form of r] corresponding 
to these three cases in i|5.4.1l to 8^5.4.31 and combine them in 
31X4] 



By approximating equation (|53p with Jg^e + J,;,e = 0, 
we have 



Q, = - 



eUgC 



where 1 + x = 1 + 



■Kn^rs^ (1 + x) 



(112) 
(113) 



5.4-1 Analytic formulae for charge-up phase / dust phase 
boundary 

We first calculate J?''^'*', the value of rj corresponding to the 
(c)charge-up phase / (d)dust phase boundary. 

The boundary satisfies J^,g = Ji,^ (Table |4)|. We use 
the approximations in (c)charge-up phase to find the break 
point of the phase. Then 



J, 



—QiU^ Tvr^'^a^^iVi 



(104) 
(105) 



because we can ignore the neutralization current j'"' , ap- 



proximate (T^ 



(geometric cross sec- 



tions) and cr^,i — -Kr^"^ \&Ql \ {"''l^bTh^)'^ (qq'a~^ » fcsT 
limit of Coulomb cross section (|31|l ). 

In dust charge-up phase, both ions and electrons are 
mainly absorbed by smaller dust, so from (|52|) and (|53|l we 
have 

_ eCug 



Qe = 



-eCn-9 



(106) 
(107) 



the factor (1 -I- x) comes from the Coulomb cross section 

Substituting and Qe, the equality for the lightning 
condition (1631) becomes 



(114) 



27rAq^,c n.Vg^ (1 + x) AWio 



rigr^ C, rUeVeU^ 



By substituting = ry^/tM^ and by solving for rj^J^, 
we have the following analytic formula for rj^'^^^: 



I 27r(l + x) Ag^„P n^VgSy^ kgT 



(115) 



We have introduced another nondimensional correction con- 
stant a''*' as we did in ^5.4.11 

5.4-3 Analytic formula for ricvit in charge-up phase 

Finally, we derive the analytic formula of the critical density 
'?criu where the condition for electric discharge (|63|l is met 
in (c)charge-up phase (c.f. H5.2.3I Figure (S^c)). 

By approximating equations (|52|l and (|54p with Js,i ~ 
Ji^e and Q„ + Qi = 0, we have 



and the absorption cross sections are geometric 



J. 



By substituting ([T06ll into (fTOSl 

— cC^Tlg — Je.i 



(108) 



cf. Fig. EJc) and in equation (|84|l . 

We have come to a simple result, that r^crit satisfies 



(109) 

Substituting equations (|104|l , (|108|) . together with = 
^(cd)^><i ^ ri^'^'^^n^ into equation (I109P and solving 

for dust number density r), we have 



A(7^,c n'^n'^rsr^ Av^^^ 



(110) 



We have introduced a nondimensional correction factor 
a^^''^ , a constant that does not depend on , r^^ , , Dj^ . 
We need this to compensate the error arising from using the 
formulae in (c)charge-up phase to find the break point of 
itself. The actual value for a''"^' is in ^5.4.41 



5.4.S Analytic formula for ricrit in dust phase 

Next, we derive the analytic formula of the critical density 
'7crit' where the condition for electric discharge (|63|l is met 
in (d)dust phase (c.f. Figureigd)). 

Imposing J^^g =0 in (|44|l . and by approximating the 
charge neutrality (|54|) with + = 0, we have 



2r„ 



(111) 



(116) 



In equation (|50p . we can ignore J^,e and further ignoring 
the second term in equation (|44p . we have 



(117) 
(118) 



2r 

— ^Ag^.^rign^ Aw^^svrr^^ - Qin^a^AV^ = 



where a, i = 



here we used the qq a~ » fcsT limit of Coulomb cross 
section (|3ip . 

Solving this for , we have 



And from equation (|53p we have 
_ engC 

Ve — 5 

By substituting these and Qe to the equality for the 
lightning condition (I63|l . replacing = 'y'^t^s' ^'^id = 
virit^'^j and by solving for 77^^'^, we obtain the following 
analytic formula for rj^%: 



(119) 



(120) 



'7crit 



27r Ag^,^ n^3n^r/ Aw^.^u^ (fc^T)' 



(121) 



We have introduced a third nondimensional correction con- 
stant a''^-' as we did in previous sections. 
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5.4.^ The combined analytic formula for j^crit 

The critical number density r^crit is either of virltyV'^''''\virif 
To choose the correct one, we have to consider the phase 
boundary conditions (c.f. Table |4]). Instead, we propose the 
following convenient scheme to choose the correct one: 

(d) . „ (d) ^ (cd) 

(cd) ■ r (c) ^ (cd) ^ (d) 

— ril'^l^ otherwise. (122) 

This scheme is based on the intuition that the (cd)phase 
boundary is included in both (c)charge-up phase and (d)dust 
phase. We can argue that if rj^.'^^^. > r?'-'^'^-', the number den- 
sity ry^'^'^' is not large enough to cause lightning, and that if 
virit < V^'^'''\ the number density is already large enough 
to cause lightning. 



Now, without the correction, e.g. with a 



1, the analytic values for ?7crit differs from the nu- 
merical values J)^";™' calculated in ^5.31 because of approx- 
imations used. For example, substituting the reference pa- 
rameter of Fig. m = 1.0 X 10"^ cm, = 1.0 x 10^ cm, 
-Dg — 2.0, and = 2.368, equation p22p evaluates to 

T/crit 



1.47 X 10 . For the reference parameter of Fig. [T] 



rg = 1.0 X 10" 



cm, — 1.0 cm, = 2.665, and 



= 3.0, equation (|122p evaluates to r^crit = 9.99 x 10^. 
The results of the simulations for those two parameter are 
viTr'' = 7.37 X 10^ and t?^;';™' = 3.01 x 10^ respectively. 
The analytic and simulational values agree upto a factor of 
three. 

We set the values for q^''', a''"'\ a*"*' by the condition 
that the following squared-error integral over the entire pa- 
rameter ranges (|97ll03p is minimum: 



dr^ dr^ dD^ dD^ log^a Vcrit - M 



(rtum) 
5 10 '7crit 



Jed) 



= 3.3 X 10" 



(123) 



This gives a'-"'' = 9.4 x 10"\ 
8.5 X 10"^. Taking these corrections into account, the values 
for rj'-''\ri'-'"^\ri'-'^^ are as follows: 



1.1 X 10'' 



( Ma,C ( 

V6.2 X 10^ ej \ 



n„ 



8.8 X 10"^ cm- 



4.7 X 10^=* cm-3 
4.0 X lO-^"* cm-3 



(t 
(r 



X 10"'' cm 

c 



X 10-1® sec-1 
3.4 X 10^ cm sec^i 



( 

U5.4eVy U.7 



3.4 X 10^ cm sec-i 



,(cd) 



3.3 X 10^ 



X 10^ K 



(124) 



6.2 X 10^ e 



n„ 



8.8 X 10-1 cm-3 



4.7 X cm-3 



4.0 X 10-1'' cm-3 



.1.0 X 10 



cm/ V 1.0 cm/ 



c 



1.0 X 10-1** sec- 
3.4 X 10^ cm sec-i 



(125) 



(d) 



5.9 X 10 



*( 



6.2 X 10^ e 

-1 



( 



4.7 X 10" cm- 



8.8 X 10-1 cm-3^ 

( Is ( \ 

V 1.0 X 10-'' cm/ V 1.0 cm/ 

c 



( 



1.0 X 10" 
15.4 eV 



3.4 X 10^ cm sec-i 



(it 



X 10^ K 



(126) 



Note that Ag^.^, n^, and also depends on dust 
parameters: , r^, , and . Using equation (|38p and 
the rich, Och introduced in i)3.4.1l 



6.2 X 10% ■ ^"'^ 



O'ch 



0.1 6.2 X 10^ e cm-2 



1.0 X 10-^ cm2 
Using equations (I5I11|I and equation ([2^ 



(127) 



(128) 



4.0 X 10' 



(-)■ 



(2.7 Au) 



-11/4 



.2.7 AU 
nT = 4.0 X 10 



\rmJ 



(2.7 Au) 



-11/4 



3.9 X 10- 



.2.7 AUy V 3.9x10-1^ 
and the monomer radius 
= 1.0 X 10-^ cm. 



(129) 



(130) 



(131) 



We have plotted these analytic solutions (|124ll26p com- 
bined with the condition (|122p in solid-line contours from 
Fig.[6]to Fig.[8l The red thin contour represents the param- 
eter ranges where r;''^''^ contributes. The blue thick contours 
represents the parameter ranges where r;'''' contributes, 
where blue solid contour means x < 1 s-nd blue dashed con- 
tour X > 1- The thick yellow-sleeved red contours represents 
the parameter ranges where rf-'^'' contributes. The numerical 
solutions, on the other hand, are plotted in colour maps and 
the black dashed contours. 

The averaged plots. Fig. [S] shows the agreement of 
the numerical and analytic value over the entire parameter 
range. Quantitatively, the root-mean-square error is 



//// dr, dr^ dD, dD^ ( 

log 10 '7crit — log 10 Vcrit 



JJJJ dr, dr^ dD, dD^ 



= 9.2 X 10" 



(132) 

Moreover, using the values of obtained only 

from the 'normal' run (Fig. |6] [7] and [8}, we can fit the 



16 Takayuki Muranushi 



results of the 'pessimistic' simulations (Fig. |9} by a root- 
mean-square error of 2.6 x 10~^. We also perform the sim- 
ulations with smaller values of relative velocity and fit the 
results. The root-mean-square errors were 5.6 x 10~^,6.4 x 
10"^l.l X 10"\ for Au^,s =u^= 3.4 x 10^ cm s-\3.4 x 
10^ cms~^,3.4 cm s~^, respectively. These fits prove the 
predictability of our analytic formulae (|122p and (|124ll26p . 



6 CONCLUSIONS AND DISCUSSIONS 

We have shown that as dust number density 77 increase, the 
charge density distribution experience four phases: (a)ion- 
electron plasma phase, (b)ion-dust plasma phase, (c)charge- 
up phase and (d)d ust phase. The former two phases are 
studied in detail by lOkuzum 1 (|2009l ). while the latter two 
phases are unique results of taking dust-dust collision into 
consideration. We have calculated the dust number density 
?7crit at which lightning strikes, as function of dust radius 
, and fractal dimension , numerically. Using the 
numerical results we have derived the analytical formulae 
for r^crit: equations (|122|l . (I124ll26p . Because the generated 
ele ctrostatic field - Kmax( ?7) grows more rapidly than estimate 
bv lGibbard et all (|l997h in (c)charge-up phase and (d)dust 
phase, lightning in protoplanetary discs are possible with 
smaller dust number densities. We discuss the consequences 
in this section. 



6.1 Energetics and direct observations 

We estimate the total energy of a lightning event in a pro- 
toplanetary disc at r = 2.7 AU. For MMSN, the num- 
ber density of the gas is 4.7 x lO'^^ cm"'' in the region. 
The typical electron mean free path at this site is Imfp — 

1.2 x lO'^ cm. By equation (|56p we know the critical elec- 



tric field Edis — 4.3 x 10~ G. The sphere with radius of 
the disc scale-height h ~ 2.4 x lO^'^ cm contains the electric 
energy W = E^-J x 47r/i^/3 ~ 4.3 x 10^^ erg. When the 
lightning strikes, the energy is concentrated into lightning 
bolt of radius w and length h, where w is related to Imfp by 
w ~ 5000 Imfp ^ 6.0 X 10^ cm (iPilipp et al.lll992l ). If all the 
energy is used to heat the gas within the lightning bolt, the 
gas can be heated to 1.6 x 10^ K. 

The ultimate energy source for this electric discharge 
event is the gravitational energy of the accreting matter. In 
our model the mass accretion ratio of uncondensed larger 
dust is M = 27rrE^u^ ~ 3.3 x lO" g sec"\ The grav- 
itational energy released within condensation region h is 
L = GMQMhr~^ ~ 6.6 x 10^* erg sec~\ For the largest 
energy event W = 4.3 x 10^^ erg. The upper limit of the 
event rate is 1.5 x 10~^ sec~^. 



modelled as solutions of Maxwell eq uations, including light - 
ning current as a source term (e.g. iRakov fc Umanlll998l ). 
When we apply these models to the protoplanetary discs, 
the electromagnetic wave spectrum is extend between the 
event duration and light crossing time of the system, or 
9.6 X 10~^ ~ 1.2 X 10~^ Hz. This frequency range is at least 
two orders of magnitude lower than any frequencies with 
established observational methods. It is difflcult to make a 
fair choice for the successor to the frequency list 'very low 
frequency (VLF),' 'ultra low frequency ( /7Z/F),' 'super low 
frequency {SLF),' and 'extremely low frequency (ELF).^ We 
opt for Astronomically Low Frequency (ALF) waves and 
hope that the reader will forgive us! Anyway the frequency 
is so low that we will need an astronomical budget to build 
an astronomically large detector to receive it, considering its 
wavelength of order of an astronomical unit. 



6.1.2 Infrared (IR) Observations 

The energy of the lightning contributes to the local heating 
of the protoplanetary discs, which might be resolved by ad- 
vanced telescopes such as Atacama Large Millimetre Array 
(ALMA). The most possible observational evidence is excess 
of heating near the snowline. To distinguish the cause of the 
heating with other heating model candidates, the variabil- 
ity or correlation function of the heating might be useful. 
This is because lightning propagates at the speed of ionised 
electrons, which is much faster than the speed of sound. 



6.1.3 Ultraviolet (UV) Observations 

The ionisation electrons of the lightning excite various elec- 
tron levels in gas molecules and dust. There is possibil- 
ity of observing fluorescence photons from such excited 
molecules. Although the disc gas is generally expected to be 
thick for ultraviolet photons, there are categories of light- 
ning that extends towar d thin regions o f the gas, known 
as sprites and elves (e.g. IWilliamg |2001^ . The sprites and 
elves are phenomena similar to lightning observed in the 
mesosphere of the earth, possibly caused by electric fields 
induced by the thunderclouds. Fluorescence lines from such 
regions c an be observed by future ultraviolet missions like 
THEIA (jSpergel et al.|[2009l ). Also, some observational re- 
sults on protostellar and protoplanetary systems today have 
diffi culties in e xplaining either lack or excess of UV (e.g. 
Nomura fc Millar 2005: iChapiUon et al.1 12OO8I : iPerez et al.l 
20081 : iHerczeg fc Hillenbrandll2008D . If excess of i/F photons 



is observed compared to the model, it might be from the 
sprite discharges and elves from the surface of the protoplan- 
etary discs; on the other hand if the chemical composition 
model require more f/F photons than is observed, lightning 
hidden in the disc mid-plane might be providing them. 



6.1.1 Astronomically Low Frequency (ALF) Waves 

The change density evolution, electromagnetic pulse, and 
electromagnetic waves accompanying li ghtning in terrestrial 
thunderclouds are observed (e.g. iKoshak fc Kriderl 1 19891 : 
iLin et al.|[l979l ). The typical wavelength of the electromag- 
netic waves are similar to the scale height of the thun- 
dercloud. These are called extremely low frequency waves. 
The electromagnetic waves from lightning can be basically 



6.1.4 High Energy Gamma Rays 

Detection of burst-like gamma-ray is reported from terres- 
trial thunder clouds. The burst precedes a cloud-to-ground 
lightning, lasts for ~ 40 seconds, extends to 10 MeV. The 
spectrum can be interpreted as consis ting of bremsstrahlung 
photons from re lativistic electrons (|Tsuchiva et al.l I2OO7I : 
lEnoto et al.ll200^ ) . These relativistic electrons are secondary 
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electrons generated by cosmic rays, and accelerated by the 
electr ic fields througli process known as avalanche amplifica- 
tion (jRoussel-Dupre fc Gurevichl *1996'). If a charged parti- 
cle is accelerated by the protoplarietary thundercloud field, 
through similar process, its kinetic energy reaches eEh ~ 

3.1 X 10^^ eV. 

6.2 Chondrule heating by lightning 

Chondrule hea ting by lightning sc e nario is now consid- 
ered unh kely ^ Weidnschillin3 1 19971 : iGibbard et al] 1 19971 : 
iGiittler et al. 200^ 1. The reasons that prohibit the scenario 
can be summarized as following three problems. 

6.2.1 Energetics problem 

The ultimate energy source (gravitational potential of the 
protoplanetary disc) is sufficient to melt the chondrules; 
but most of the energy earned by ingoing larger dust 
go to the outgoing gas by angular momentum exchange 
l| Weidnschillind [l997l ) : little contribute to the random mo- 
tion, the energy source for the lightning. 

6.2.2 Neutralization problem 

Unlike the earth atmosphere, the protoplanetary discs are 
filled with weakly ionised plasma which rapidly responds to 
electric field. Neutralization effect can be further subdivided 
to microscopic neutralization of individual dust and macro- 
scopic neutralization of large-scale electric field necessary to 
cause lightning. If a dust get charged by dust-dust collision, 
the dust absorbs plasma of opposite polarity in ~ 10 sec 
and returns to equilibrium charge state. Moreover, even if 
there is charged dust and bulk motion between the oppo- 
sitely charged dust, the electric field caused by the dust in- 
duces Ohmic current in the plasma. The current will quickly 
neutralize the electric field. 

6.2.3 Destruction problem 

After all, there is an experimental evidence bv lGiittler et al.l 
l|2008t ) that lightning destroys the dust aggregates rather 
than melting them. 

6.2.4 Solution to the problems 

This work can provide answer for the first and second prob- 
lem, energetics problem, the larger dust and the gas (con- 
taining smaller dust that are coupled to the gas) is now 
'harnessed' by electric field. Outgoing gas is not free in car- 
rying the gravitational energy away; instead the gas con- 
verts its gravitational energy into electric field energy, fully 
contributing to lightning. For the neutralization problem, 
we have shown in this work that with reasonably high dust 
number density rj, the dust-dust charge separation can dom- 
inate over the plasma neutralization, and the electrostatic 
field can grow up to critical value. 

For the third problem, we point out that in 
IGiittler et al.l (|2008l )'s experiment, either the electron mean 
free path is many orders of magnitude shorter, or the 
electron kinetic energy is much larger compared to the 



protoplanetary-disc environment. They used air at pressures 
between 10 and 10^ Pa. Air consists of 78 per cent nitrogen, 
21 per cent oxygen, and 1 per cent argon. Their molecular 
van der Waals radii are 1.6 x 10~^ cm, 1 .5 x 10~* cm, and 
1.9 X 10"** cm, respectively (|Bondilll964l ). 

Therefore, the electron mean free path and the electron 
kinetic energy, We = eSZ^fp, was Zmtp ~ 4.8 x 10~^ cm, 
We = 1.6 X 10" eV for 10 Pa case, and Z^fp ~ 4.8 x lO"'^ cm , 
We = 1.6 eV for 10'' Pa case, respectively. On the other hand 
in protoplanetary discs, typical mean free path and electron 
kinetic energy are In-iip = 1.2 x lO'^ cm. We = 15.4 eV. 

It might be possible that protoplanetary-disc lightning 
is effective in melting dust aggregates, although experimen- 
tal lightning is ineffective in heating and led to disruption 
of the dust, due to shorter mean free path or higher energy 
electron. The minimum size of the structures that electron 
can form is of order of its mean free path. If the electron 
mean free path is much shorter than the dust aggregates, 
as in 10^ Pa case, the electron current may concentrate on 
the most conductive part of the dust aggregate, leading to 
partial heating and explosion of the dust. On the other hand 
if the electron is much more energetic, as in 10 Pa case, it 
may react differently on dust monomers. 

To reproduce the mean free path and electron energy si- 
multaneously, one must reproduce the electric field strength 
E = 4.3 x 10"** G of protoplanetary discs; while the elec- 
tric field used in the experiment i5 = 1.1 x lO'^ G was 
much stronger. This much stronger electric field itself, might 
be the cause of dust aggregate dissociation, due to much 
stronger electric force exerted on electron-absorbed dust 
monomers. Also the discharge time-scale in the experiment 
was much smaller than that in the protoplanetary discs, 
which might have led to the catastrophic results. 

We think that the effect of lighting on dust aggregate 
in protoplanetary-disc environment is yet to be confirmed in 
future experiments and simulations. 



6.3 Effects on magnetorotational instabihty 
(MRI) and disc environment 

The dust-dust collisional charging and lightning is not a side- 
effect of some other processes, but is one of the key processes 
in protoplanetary discs that affects each other. The lightning 
is powered by gravitational energy of the migrating larger 
dust. The migration of the larger dust as well as the long 
term evolution of the gas disc is governed by the disc vis- 
cosity. The best candidate for providing the disc viscosity is 
MRI. And MRI is controlled by gas ionisation degree, which 
in turn is controlled by the dust charge state and lightning. 

Even the longest estimate for time-scale of the light- 
ning 1.0 X 10* sec is much smaller than the time-scale of 
MRI, which is at least of the order of Kepler timescales. 
Lightning occur in low-ionisation regions where MRI is pro- 
hibited (dead zones), and even if the lightning instantly raise 
the ionisation rate, the free electrons and ions will quickly 
be absorbed by the dust. Therefore we expect that MRI and 
lightning cannot co-exist. However lot of profound phenom- 
ena are possible. Just for an example let us think of a two- 
layer dead-active zone model but with dust-dust collisional 
charging. The dead-zone is filled with lightning, inducing 
sprite discharges towards active zones, which sustains the 
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ionisation rate and MRI. The MRI in turn shovels the dust 
into the dead-zone. 

Such global models are beyond the reach of this pa- 
per. Nevertheless we conclude this paper by stating that 
the dust-dust coUisional charging is a necessary component 
for understanding the planetesimal formation and global be- 
haviour of the protoplanetary discs. 
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APPENDIX A: CATIONIC DUST AND 
ANIONIC DUST 

In this section we justify the two-dust picture introduced 
in >i2.1l Protoplanetary discs consist of dust with various 
parameter J. The parameter vector J may include, but is 
not limited to, dust radius, porosity, material and surface 
chemical potential. We classify these dust into two groups 
according to their electric tendency; one is cationic dust who 
receive positive charge through dust-dust collision and the 
other is anionic dust who receive negative charge. In this 
section we give the precise definition of cationic and anionic 
dust. 

Let rij indicate the number density of the dust with 
parameter J. Let (Tj'^j, Avy^j, Aqy j be collisional cross sec- 
tion, mean relative velocity, and mean amount of charge 
that moves from dust J' to dust J in a collision, respec- 
tively. Then we can calculate qj, charge received by dust J 
per unit time, by 



j'j"'j'<jj'j 



(Al) 



We define cationic dust and anionic dust as C = {J\qj > 
0} and A = {J\qj < 0}. Cationic dust receive net positive 
charge in dust-dust collision and tend to be cationic, while 
anionic dust tend to be anionic. 

We assume the average dust distribution rij"' as that of 
MMSN model. We also assume that at some local condensa- 
tion regions, dust number density is multiplied by carrying 
in dust from other portions of the disk. For simplicity we 
assume that the relative number density rj is independent 
of dust parameter j so that nj = ryrij"' (e.g. this is the case 
when collisional cascade equilibrium is faster than the mi- 
gration). We define ayj and Aq/ j as the value for neutral 
dust. As the dust acquires charge, the amount of charge 
exchanged in a single collision becomes smaller due to the 
exchange of the charge they already have; we treat this devi- 
ation from neutral dust as separate 'neutralization' channel. 
With these assumptions, the sign of qj HA1|) do not depend 
on rj, and the term 'cationic dust' and 'anionic dust' is well 
defined independent of dust number density ry. 

Now we can simplify the problem by treat cationic and 
anionic dust as if they are two discrete kinds of dust. There- 
fore we define the representative variables for cationic and 
anionic dust as follows. 



nj = rj n[°^ 



E 



JSC 

n'j 

J6A 



jjeA,j'ec 



njnjiVjf 



njUjiVj^ji 



jeA,j' gc 



Aqj^j>aj^j, UjUji Vjy 



(A2) 
(A3) 

(A4) 
(A5) 
(A6) 
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Figure 6. Value of fjcrit as function of dust radius r^, and fractal dimension D^, Dj^. The base values are = 1.0 X 10~^ cm, 
= 1.0 X 10^ cm, rrig = 3.9 X 10~^ g, = 1.5 X 10^ g, Dg = 2.0, and = 2.368. We keep mi constant when we vary Di\ we 
keep Di constant when we vary ri. Numerical results are in colour maps and black dashed contours; analytical values in coloured solid 
contours (c.f. ^5.4.4! for the details of the plots.) 
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The base values are r„ 



1.0 X 10' 



= 1.0 cm, rrig = 1.9 X 10^^^ g, rrij^ = 3.9 g, Dg = 2.665, and = 3.0. We keep mi constant when we vary Di\ we keep Di constant 
when we vary Ti. Numerical results are in colour maps and black dashed contours; analytical values in coloured solid contours (c.f. N5.4.4I 
for the details of the plots.) 
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Figure 8. Value of r^^rit ^ function of r^, r^, , and . Parameters do not appear in x-axis or y-axis arc uniformly averaged over the 
parameter range we accept. Numerical results are in colour maps and black dashed contours; analytical values in coloured solid contours 
(c.f. i|5.4.4l for the details of the plots.) 
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Figure 9. In the above plots, all the parameters but the charge separation efficiency is same as that of Fig. |6] while the the charge 
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Dust-dust Collisional Charging and Lightning 25 



APPENDIX B: SIMULATIONS 

In this section, we briefly describe our numerical methods. 
We need to solve the equilibrium equations (|50ll54p , for vari- 
ous environmental parameters. Especially we vary r) for each 
set of other parameters. Then we know the minimum rj that 
satisfies the electric discharge condition (|63|l . for the each 
set of other parameters. 

This kind of problem, a massive parameter parallelism, 
is typicall y suitable f or massively parallel computing hard- 
ware (e.g. lFord|[2009l ), such as gener al purpose gra phic pro- 
cessors (GPGPUs) or GRAPE-DR l|Makindl2008l ). We de- 
scribe the CPU and GPGPU based programmes we used in 
this research to solve equations (|40I43|I in this section. 



Bl Direct integral solver 

The most straightforward means of finding the equilibrium 
solutions (|4Uj|43p is to directly integrate the dynamic equa- 
tions (|40l43p . Nevertheless, constant-time-step direct inte- 
gral cannot solve (|40I43[) correctly for some of the parameter 
range. This is because the current densities Qi differ many 
orders of magnitude for such parameters. We must choose 
the integration time-step dti carefully. This leads us to the 
use of a adaptive time step. 

The adaptive-time-step direct-integral solver follows the 
dynamic equations (|40I43|) in terms of discretized time ti 
where time ti is incremented by dynamic time-step dti: 



Qi,i+i — Qi,i + Ji' ,ldti, 
i' 

ti+i = ti dti. 
Our choice of the dynamical time-step dti is as follows; 



ratio (I) 



dti 



Ql.i ~ — 1 



1.0 X 10"'' • dti^i 
maxi {ratio (I)) 



(Bl) 
(B2) 

(B3) 
(B4) 



The direct integral solver is reliable, in sense that it 
is less prone to implementation mistakes because it almost 
straightforwardly reflects the equations (|40l43p . and that out 
of possible many equilibrium solutions (|40I43|I the solver will 
always find the desired equilibrium. 

However, as rj become much larger or much smaller than 
unity, we have found that charge distribution in the system 
get unbalanced. As we try to update the species I with lit- 
tle charge but large current, the dynamic time-step dti (|B4|) 
become the time-scale the equilibrium is reached, and the 
simulations becomes time consuming. Use of higher-order 
integral schemes are futile because we cannot take time- 
step much larger than dti. In addition to that, floating point 
numbers mainly available on GPU are single-precision float- 
ing point numbers. Computation of double-precision floating 
point numbers are either not supported or order of magni- 
tude slower on common GPU. 



B2 Binary search solver 

So, we need to flnd an alternative method to solve the equi- 
librium equations (|50II54|) without directly integrating the 



dynamic equation, avoiding the addition between numbers 
of different magnitude as long as possible. 

Binary search is a common method for flnding zero 
point of a one-parameter function /; solving equation f{x) = 
for X. To solve the system of equations (|50j|54|l . we divide 
the problem into set of one-parameter problems, and con- 
quer by recursive use of binary search. In doing so, we have 
to be careful in choosing which of equations (|50j|54p we use 
to flnd zero point of which freedom Qi. Wrong choice leads 
to wrong result. 

First, Qi and Qe can be analytically expressed in terms 
of and Qj^ as follows: 



Where we made abbreviations 

f^L.i = '^cou{qL,e,r^,kBT) , 
i^L.e = (ycou{qL,~e,r^,kBT) 



(B5) 
(B6) 



(B7) 
(B8) 



and so on. Further eliminations of freedoms is possible but 
complicated, because of complicated and sign-sensitive form 
of the Coulomb cross sections (13013111 . Instead we resort to 
numerical methods to flnd out the equilibrium solution 
and Qj^ for each rj, and then flnd T^crit, the minimum rj that 
satisfles the electric discharge condition (|63p . 

We now describe how to solve the system of equations 
(|50p . (|5ip . (|54p . and to flnd ?7crit, provided that for any one- 
parameter / we can solve f{x) — 0. 

Let us name the left-hand-sides of equations (|50p . (|5ip . 
(|54p as fj^, /s , and /s . We eliminate Qi and Qe from 
these using (|B5P and (|B6p . and regard them as functions 
of 

'?! Ql ' Qs ^ follows: 



fs {v,Ql^Qs ) 

/crit {V,Ql'Qs) 



-J, 



J. 



J. 



^L,S ^S^' "'S.e' 
= Ql +Qs +Q^{VlQL^Qs) 

+ Q^{v,Ql'Qs)' 

_ Ql rrieVeUL 

QeiV,QL,Qs) A Won' 



(B9) 
(BIO) 

(BU) 
(B12) 



We have also deflned fait according to (|63p . 

For each flxed set of rj and Qj^ , /s Q^ , Qs ) is a one- 
parameter function of Q^ . According to the assumption we 

can solve /e (??, Q^, , Qs ) ~ fo'" Qs ■ We deflne (v^Ql) to 
denote the solution, so that 



h [V, 



holds. 



(B13) 



Then for each flxed rj, f^irjjQ^TQ'^irjjQL)) a one- 
parameter function of Qj^ . According to the assumption we 
can solve {rj,Q^,Ql{rj,Q^)) =0 for Q^ . We deflne Q° (rj) 
to denote the solution, so that 



/. {v,Ql {v),Ql {v,Ql {ri)))=0 



(B14) 



holds. 



Then /c,it('), Q^J (??), Q° (??))) is a one-parameter 

function of rj. According to the assumption we can solve 
/crit(?7,Q2(??),Q°(r?,Q°(77))) = for rj. We deflne rj° to de- 



26 Takayuki Muranushi 



note the solution, so that 

L{n\Ql{n'),QlW,Ql{ri")))=o (bi5) 

holds, which is the r^crit we are looking for. 

With this method, whenever a solver approaches the 
zero point of one of /, the / will consist of two or more terms 
of same magnitude, and of other terms of smaller magnitude. 
Smaller terms are irrelevant to the equilibrium. So we will 
always be comparing the terms of same magnitude. Because 
of this, the method gives sufficiently precise solutions even 
with single precision floating point numbers. 



files and link the object files using NVCC, CUDA compiler pro- 
vided by NVIDIA. 

Another example of such benefit is thrust 
(http://code.google.eom/p/thrust/), a CUDA coun- 
terpart of what standard template library (STL) is in 
CH — h. With THRUST, for example, device and host memory 
management is automated. Memories are allocated and 
freed automatically in the constructor and destructor of 
the container classes. Copying data between host memory 
and device memory are simply expressed by substitution = 
operators. 



B3 Binary search solver on GPGPU 

Graphic processing units (GPUs) are processors special- 
ized for computer graphics, widely used in personal com- 
puters, workstations, and video game devices. But as more 
and more realistic computer graphics have been demanded, 
GPUs became capable of more and more types of calcu- 
lations, and finally evolved into general purpose graphic 
processing units (GPGPUs), who are programmable for 
general computation, not limited to graphic processing. 
Due to the nature of graphic processing tasks, GPUs' 
parallel computation capacities are are one or two or- 
ders of magnitude larger compared to that of CPUs. On 
the other hand marketplace competition and mass produc- 
tion keep the CPUs' price low. Although parallel program- 
ming has been a hard task for programmers, the paral- 
lelism found in nature, together with GPGPU's power and 
price makes it very alluring as next generation computa- 
tional platform for computational astrophysics, and com- 
putational natural science. Use of GPGPU have already 
started in several fields of astronomy and astrophy sics, such 
as si gnal processing (e.g. iHarris et al . 2008: Wavth et al. 
200£ ) , N-body simulations of gravitv (e.g. [Hamada fc litaka 
2007l : lBeUeman et al.ll2008l:lMoore fc Qunienll2008l '). gravita- 
tion al lensing ( e.g. Thompson et al.l 2009|). orbital dynamics 
(e.g.|Fordll2009 ^). radiation-transfer (e.g. Jonsson fc Primackl 
20091'). and al so in var ious other bran c hes of science (e.g . 
van Meel et a l. 2007; AndrecutI l2008l : iBarros et all l20oi 
Januszewski fc Kostuill2009l ). 

With GPGPU we can challenge problems that had 
been computationally formidable. To begin this challenge, 
we have constructed Tengu, (Tenmon-GPGPU cluster; 
GPGPU cluster for astrophysical purposes. It consists of 10 
computer nodes, each equipped with NVIDIA's GPGPU. 
We use the programming language CUDA to write codes for 
GPGPUs. CUDA is compatible with CH — h, so we benefit both 
from expressive power of CH — h and computational power of 
GPGPUs. 

Thanks to this, we organize our CH — h and CUDA codes 
in the following way. We made CH — h classes representing the 
protoplanetary disc, dust plasma, problem initial conditions 
and solutions, and the numerical solvers. Each solver inherit 
from an abstract solver class. Thus the users of the solvers, 
namely the programme parts that carries out tests and nu- 
merical experiments can use any solver they prefer, without 
noticing what algorithm the solvers use nor on what hard- 
ware they run. We write most of the code in CH — h and com- 
pile them by GCC. The GPGPU related details are separated 
in several .cu files by means of pimpl idiom. We compile .cu 



B4 Testing 

We choose the Test-Driven Development style for this study. 
We have tested that the charge conservation and the current 
conservation conditions are held, for each equilibrium solu- 
tion that each solver give. We have also tested that the value 
of the currents satisfy equations by comparing them with the 
simplest implementation. 

Why do we test our codes? We need tests in numerical 
physics because the codes must compile correctly, the codes 
must translate the algorithms correctly, and algorithms rep- 
resent the physical concepts correctly. 

In order to assure these, we are accustomed to perform 
various tests during a code development, by examining the 
internal states and outputs of the programme. Furthermore, 
we want to make sure that criteria once tested always meet 
thereafter, and that the tests cover all the important aspects 
of the code. As the code grow, it becomes more and more 
effective to build up a system of test rather than to perform 
tests manually. This technique is know n as Test-Driven De- 
velopment (e.g. lErdogmus et al.ll2005l '). 

The programme is divided into many functional mod- 
ules, and we test that these modules give expected out- 
put for given input. This is compared to code-reading 
style of tests, where testers finds faults in the code by 
reading them. Alth ough code reading is capable of find- 
ing more mistakes (|Basili fc Selbvl 1 19871 ) it is only effec- 
tive when the code is short and it is possible to trace 
the comprehensive behaviour of the code line-by-line. Unit 
tests, on the other hand cares only on the input and out- 
put. It is effective even if we are trying new languages 
or hardware, or we cannot debug trace on them. We use 
GOOGLETEST, Google's framework for writing automated 
CH — h tests | |http : //code .google . com/p/googletest7] ). 

We must also consider the time cost of the test. Because 
we do computationally heavy tasks, examining the entire 
behaviour of the programme is not practical. With system- 
atized tests, we can ensure the equivalence of the codes as 
we optimize them, or as we transplant it onto another lan- 
guage or hardware. We further construct a 'test ladder,' an 
analog of distance ladder in cosmology. We develop a chain 
of successively faster algorithms, and feed them with ran- 
domly generated inputs and check if their response is same 
up to required precision. At the one end of the ladder is a 
code that almost directly trace the equations, correctness of 
whose implementation is self-evident. At the other end of 
the ladder is optimized, massively parallel code running on 
GPGPU. 

Finally, we evaluate the computational optimization 
achieved by measuring the speed of each solvers in terms 
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solver 


# of cases 


time 


optimization 


direct integral, CPU^ 


9885 


636714 


1.0 


binary search, CPU^ 


19200 


977.973 


1.3 X 10^ 


binary search, GPU^ 


19200 


6.98705 


1.8 X 10^ 


final problem* 


364325 


226.469 


1.0 X 10^ 



Table Bl. The name of the computation runs, the size of initial 
conditions sets, and the wall clock time needed to solve »7crit for 
all initial conditions in seconds. The speed of the codes are also 
listed, in terms of number of solved cases per time, relative to the 
first case. (1) We used the direct integral solver (||BTJ, running it 
in parallel on 40 CPU cores. (2) We used the binary search solver 
(gHU, on a Core 2 Quad 9300 CPU in single thread. (3) We used 
the CUDA version of binary search solver, ( i|B3l l. on single GTX280 
GPGPU. (4) We used the same programme ( i|B3l l and the same 
GPU for actual numerical experiment. The run generates a set of 
data that corresponds to rjcrit as function of dust radius , 
and fractal dimension , . Or it corresponds to one page of 
the result figure, e.g. Fig. |6] 

of how many problems they solve per wall clock time. See 
Table iBll for optimization results. 



